Load all required libraries.

library(tidyverse)
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v ggplot2 3.3.5     v purrr   0.3.4
## v tibble  3.1.3     v dplyr   1.0.7
## v tidyr   1.1.3     v stringr 1.4.0
## v readr   2.0.0     v forcats 0.5.1
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(broom)

Read in raw data from RDS.

raw_data <- readRDS("./n1_n2_cleaned_cases.rds")

Make a few small modifications to names and data for visualizations.

final_data <- raw_data %>% mutate(log_copy_per_L = log10(mean_copy_num_L)) %>%
  rename(Facility = wrf) %>%
  mutate(Facility = recode(Facility, 
                           "NO" = "WRF A",
                           "MI" = "WRF B",
                           "CC" = "WRF C"))

Seperate the data by gene target to ease layering in the final plot

#make three data layers
only_positives <<- subset(final_data, (!is.na(final_data$Facility)))
only_n1 <- subset(only_positives, target == "N1")
only_n2 <- subset(only_positives, target == "N2")
only_background <<-final_data %>% 
  select(c(date, cases_cum_clarke, new_cases_clarke, X7_day_ave_clarke)) %>%
  group_by(date) %>% summarise_if(is.numeric, mean)

#specify fun colors
background_color <- "#7570B3"
seven_day_ave_color <- "#E6AB02"
marker_colors <- c("N1" = '#1B9E77',"N2" ='#D95F02')
#remove facilty C for now
#only_n1 <- only_n1[!(only_n1$Facility == "WRF C"),]
#only_n2 <- only_n2[!(only_n2$Facility == "WRF C"),]

only_n1 <- only_n1[!(only_n1$Facility == "WRF A" & only_n1$date == "2020-11-02"), ]
only_n2 <- only_n2[!(only_n2$Facility == "WRF A" & only_n2$date == "2020-11-02"), ]

Build the main plot

      #first layer is the background epidemic curve
        p1 <- only_background %>%
              plotly::plot_ly() %>%
              plotly::add_trace(x = ~date, y = ~new_cases_clarke, 
                                type = "bar", 
                                hoverinfo = "text",
                                text = ~paste('</br> Date: ', date,
                                                     '</br> Daily Cases: ', new_cases_clarke),
                                alpha = 0.5,
                                name = "Daily Reported Cases",
                                color = background_color,
                                colors = background_color,
                                showlegend = FALSE) %>%
            layout(yaxis = list(title = "Clarke County Daily Cases", showline=TRUE)) %>%
            layout(legend = list(orientation = "h", x = 0.2, y = -0.3))
        
        #renders the main plot layer two as seven day moving average
        p1 <- p1 %>% plotly::add_trace(x = ~date, y = ~X7_day_ave_clarke, 
                             type = "scatter",
                             mode = "lines",
                             hoverinfo = "text",
                            text = ~paste('</br> Date: ', date,
                                                     '</br> Seven-Day Moving Average: ', X7_day_ave_clarke),
                             name = "Seven Day Moving Average Athens",
                             line = list(color = seven_day_ave_color),
                             showlegend = FALSE)
      

        
        #renders the main plot layer three as positive target hits
        
        p2 <- plotly::plot_ly() %>%
          plotly::add_trace(x = ~date, y = ~mean_copy_num_L,
                                       type = "scatter",
                                       mode = "markers",
                                       hoverinfo = "text",
                                       text = ~paste('</br> Date: ', date,
                                                     '</br> Facility: ', Facility,
                                                     '</br> Target: ', target,
                                                     '</br> Copies/L: ', round(mean_copy_num_L, digits = 2)),
                                       data = only_n1,
                                       symbol = ~Facility,
                                       marker = list(color = '#1B9E77', size = 8, opacity = 0.65),
                                       showlegend = FALSE) %>%
          plotly::add_trace(x = ~date, y = ~mean_copy_num_L,
                                       type = "scatter",
                                       mode = "markers",
                                       hoverinfo = "text",
                                       text = ~paste('</br> Date: ', date,
                                                     '</br> Facility: ', Facility,
                                                     '</br> Target: ', target,
                                                     '</br> Copies/L: ', round(mean_copy_num_L, digits = 2)),
                                       data = only_n2,
                                       symbol = ~Facility,
                                       marker = list(color = '#D95F02', size = 8, opacity = 0.65),
                                       showlegend = FALSE) %>%
            layout(yaxis = list(title = "SARS CoV-2 Copies/L", 
                                 showline = TRUE,
                                 type = "log",
                                 dtick = 1,
                                 automargin = TRUE)) %>%
            layout(legend = list(orientation = "h", x = 0.2, y = -0.3))
        
        #adds the limit of detection dashed line
        p2 <- p2 %>% plotly::add_segments(x = as.Date("2020-03-14"), 
                                          xend = ~max(date + 10), 
                                          y = 3571.429, yend = 3571.429,
                                          opacity = 0.35,
                                          line = list(color = "black", dash = "dash")) %>%
          layout(annotations = list(x = as.Date("2020-03-28"), y = 3.8, xref = "x", yref = "y", 
                                    text = "Limit of Detection", showarrow = FALSE))

        

        p1
        p2

Combine the two main plot pieces as a subplot

#seperate n1 and n2 frames by site
#n1
wrf_a_only_n1 <- subset(only_n1, Facility == "WRF A")
wrf_b_only_n1 <- subset(only_n1, Facility == "WRF B")
wrf_c_only_n1 <- subset(only_n1, Facility == "WRF C")

#n2
wrf_a_only_n2 <- subset(only_n2, Facility == "WRF A")
wrf_b_only_n2 <- subset(only_n2, Facility == "WRF B")
wrf_c_only_n2 <- subset(only_n2, Facility == "WRF C")


#rejoin the old data frames then seperate in to averages for each plant. 
wrfa_both <- full_join(wrf_a_only_n1, wrf_a_only_n2)%>%
  select(c(date, mean_total_copies)) %>%
  group_by(date) %>%
  summarize_if(is.numeric, mean) %>%
  ungroup() %>%
  mutate(log_total_copies_both = log10(mean_total_copies))
## Joining, by = c("date", "new_cases_clarke", "cases_cum_clarke", "X7_day_ave_clarke", "Facility", "collection_num", "target", "mean_copy_num_uL_rxn", "mean_copy_num_L", "sd_L", "mean_total_copies", "sd_total_copies", "log_copy_per_L")
wrfb_both <- full_join(wrf_b_only_n1, wrf_b_only_n2)%>%
  select(c(date, mean_total_copies)) %>%
  group_by(date) %>%
  summarize_if(is.numeric, mean) %>%
  ungroup() %>%
  mutate(log_total_copies_both = log10(mean_total_copies))
## Joining, by = c("date", "new_cases_clarke", "cases_cum_clarke", "X7_day_ave_clarke", "Facility", "collection_num", "target", "mean_copy_num_uL_rxn", "mean_copy_num_L", "sd_L", "mean_total_copies", "sd_total_copies", "log_copy_per_L")
wrfc_both <- full_join(wrf_c_only_n1, wrf_c_only_n2)%>%
  select(c(date, mean_total_copies)) %>%
  group_by(date) %>%
  summarize_if(is.numeric, mean) %>%
  ungroup() %>%
  mutate(log_total_copies_both = log10(mean_total_copies))
## Joining, by = c("date", "new_cases_clarke", "cases_cum_clarke", "X7_day_ave_clarke", "Facility", "collection_num", "target", "mean_copy_num_uL_rxn", "mean_copy_num_L", "sd_L", "mean_total_copies", "sd_total_copies", "log_copy_per_L")
#get max date
maxdate <- max(wrfa_both$date)
mindate <- min(wrfa_both$date)

Build loess smoothing figures figures

This makes the individual plots

#**************************************WRF A PLOT**********************************************
#add trendlines 
#extract data from geom_smooth
#both extract
# *********************************span 0.6***********************************
#*****************Must always update the n = TOTAL NUMBER OF DAYS*************************
extract_botha <- ggplot(wrfa_both, aes(x = date, y = log_total_copies_both)) + 
  stat_smooth(aes(outfit=fit_botha<<-..y..), method = "loess", color = '#1B9E77', 
              span = 0.3, n = 457)
## Warning: Ignoring unknown aesthetics: outfit
#look at the fits to align dates and total observations
#both
extract_botha
## `geom_smooth()` using formula 'y ~ x'

fit_botha
##   [1] 12.82935 12.83650 12.84355 12.85051 12.85736 12.86410 12.87071 12.87721
##   [9] 12.88357 12.88980 12.89588 12.90181 12.90759 12.91321 12.91866 12.92394
##  [17] 12.92904 12.93396 12.93868 12.94320 12.94752 12.95164 12.95554 12.95923
##  [25] 12.96273 12.96604 12.96916 12.97212 12.97490 12.97753 12.98001 12.98235
##  [33] 12.98455 12.98663 12.98859 12.99043 12.99217 12.99382 12.99538 12.99686
##  [41] 12.99827 12.99962 13.00091 13.00207 13.00303 13.00381 13.00440 13.00481
##  [49] 13.00506 13.00515 13.00508 13.00487 13.00451 13.00403 13.00342 13.00269
##  [57] 13.00185 13.00091 12.99988 12.99875 12.99755 12.99627 12.99492 12.99352
##  [65] 12.99206 12.99056 12.98902 12.98745 12.98586 12.98425 12.98264 12.98060
##  [73] 12.97776 12.97419 12.96998 12.96518 12.95987 12.95413 12.94802 12.94161
##  [81] 12.93499 12.92822 12.92137 12.91452 12.90773 12.90109 12.89466 12.88851
##  [89] 12.88272 12.87736 12.87250 12.86821 12.86349 12.85736 12.84995 12.84137
##  [97] 12.83174 12.82119 12.80983 12.79778 12.78517 12.77211 12.75872 12.74512
## [105] 12.73143 12.71777 12.70427 12.69103 12.67819 12.66585 12.65414 12.64318
## [113] 12.63310 12.62400 12.61600 12.60924 12.60383 12.59988 12.59752 12.59686
## [121] 12.59810 12.60121 12.60597 12.61217 12.61960 12.62805 12.63730 12.64716
## [129] 12.65740 12.66781 12.67819 12.68832 12.69799 12.70699 12.71510 12.72310
## [137] 12.73181 12.74115 12.75105 12.76140 12.77213 12.78315 12.79438 12.80572
## [145] 12.81709 12.82841 12.83959 12.85054 12.86118 12.87157 12.88186 12.89210
## [153] 12.90232 12.91257 12.92289 12.93333 12.94392 12.95470 12.96573 12.97704
## [161] 12.98868 13.00340 13.02335 13.04769 13.07553 13.10602 13.13829 13.17149
## [169] 13.20474 13.23718 13.26795 13.29618 13.32101 13.34158 13.35703 13.37055
## [177] 13.38575 13.40237 13.42011 13.43873 13.45794 13.47747 13.49705 13.51641
## [185] 13.53527 13.55337 13.57044 13.58619 13.60036 13.61268 13.62288 13.63068
## [193] 13.63581 13.63800 13.63698 13.63248 13.62414 13.61211 13.59686 13.57888
## [201] 13.55863 13.53660 13.51326 13.48909 13.46456 13.44016 13.41635 13.39362
## [209] 13.37244 13.35329 13.33340 13.31017 13.28432 13.25658 13.22768 13.19834
## [217] 13.16929 13.14126 13.11496 13.08734 13.05527 13.01946 12.98056 12.93928
## [225] 12.89629 12.85227 12.80790 12.76387 12.72085 12.67953 12.64059 12.60472
## [233] 12.57258 12.54073 12.50571 12.46825 12.42910 12.38897 12.34859 12.30870
## [241] 12.27001 12.23327 12.19920 12.16853 12.14199 12.11825 12.09553 12.07379
## [249] 12.05303 12.03324 12.01440 11.99650 11.97954 11.96349 11.94834 11.93409
## [257] 11.92072 11.90822 11.89658 11.88578 11.87581 11.86667 11.85833 11.85079
## [265] 11.84403 11.83804 11.83399 11.83277 11.83396 11.83715 11.84192 11.84785
## [273] 11.85453 11.86153 11.86845 11.87486 11.88034 11.88448 11.88687 11.88708
## [281] 11.89072 11.90217 11.91920 11.93955 11.96096 11.98120 11.99800 12.00913
## [289] 12.01233 12.01171 12.01263 12.01476 12.01774 12.02122 12.02487 12.02832
## [297] 12.03123 12.03326 12.03405 12.03326 12.03054 12.02554 12.01792 12.00597
## [305] 11.98885 11.96748 11.94280 11.91573 11.88720 11.85814 11.82948 11.80216
## [313] 11.77708 11.75520 11.73743 11.72018 11.69961 11.67628 11.65075 11.62357
## [321] 11.59530 11.56651 11.53774 11.50957 11.48254 11.45721 11.43416 11.41392
## [329] 11.39707 11.37953 11.35781 11.33332 11.30748 11.28170 11.25740 11.23598
## [337] 11.21888 11.20749 11.20003 11.19376 11.18868 11.18478 11.18205 11.18049
## [345] 11.18009 11.18086 11.18278 11.18585 11.19006 11.19542 11.20190 11.20952
## [353] 11.21827 11.22813 11.23911 11.25120 11.26439 11.27869 11.29408 11.31129
## [361] 11.33093 11.35278 11.37664 11.40230 11.42955 11.45819 11.48800 11.51879
## [369] 11.55034 11.58245 11.61491 11.64751 11.68005 11.71232 11.74411 11.77521
## [377] 11.80542 11.83453 11.86611 11.90293 11.94356 11.98660 12.03064 12.07425
## [385] 12.11603 12.15457 12.18844 12.22164 12.25864 12.29873 12.34122 12.38540
## [393] 12.43058 12.47607 12.52115 12.56513 12.60731 12.64699 12.68348 12.71607
## [401] 12.74406 12.76994 12.79642 12.82314 12.84973 12.87584 12.90111 12.92517
## [409] 12.94767 12.96825 12.98653 13.00217 13.01480 13.02541 13.03518 13.04403
## [417] 13.05190 13.05871 13.06440 13.06889 13.07212 13.07400 13.07448 13.07348
## [425] 13.07093 13.06676 13.06089 13.05367 13.04542 13.03606 13.02551 13.01370
## [433] 13.00053 12.98594 12.96984 12.95215 12.93295 12.91235 12.89037 12.86701
## [441] 12.84227 12.81614 12.78864 12.75976 12.72951 12.69788 12.66488 12.63051
## [449] 12.59478 12.55768 12.51922 12.47940 12.43822 12.39568 12.35179 12.30654
## [457] 12.25994
#assign fits to a vector
both_trenda <- fit_botha

#extract y min and max for each
limits_botha <- ggplot_build(extract_botha)$data
## `geom_smooth()` using formula 'y ~ x'
limits_botha <- as.data.frame(limits_botha)
both_ymina <- limits_botha$ymin
both_ymaxa <- limits_botha$ymax

#reassign dataframes (just to be safe)
work_botha <- wrfa_both

#fill in missing dates to smooth fits
work_botha <- work_botha %>% complete(date = seq(min(date), max(date), by = "1 day"))
date_vec_botha <- work_botha$date

#create a new smooth dataframe to layer
smooth_frame_botha <- data.frame(date_vec_botha, both_trenda, both_ymina, both_ymaxa)
#WRF A
#plot smooth frames
p_wrf_a <- plotly::plot_ly() %>%
  plotly::add_lines(x = ~date_vec_botha, y = ~both_trenda,
                    data = smooth_frame_botha,
                    hoverinfo = "text",
                    text = ~paste('</br> Date: ', date_vec_botha,
                                  '</br> Median Log Copies: ', round(both_trenda, digits = 2)),
                    line = list(color = '#1B9E77', size = 8, opacity = 0.65),
                    showlegend = FALSE) %>%
     layout(xaxis = list(range = c(mindate - 7, maxdate + 7))) %>% #buffer here
plotly::add_ribbons(x ~date_vec_botha, ymin = ~both_ymina, ymax = ~both_ymaxa,
                    showlegend = FALSE,
                    opacity = 0.25,
                    hoverinfo = "text",
                    text = ~paste('</br> Date: ', date_vec_botha, #leaving in case we want to change
                                  '</br> Max Log Copies: ', round(both_ymaxa, digits = 2),
                                  '</br> Min Log Copies: ', round(both_ymina, digits = 2)),
                    name = "",
                    fillcolor = '#1B9E77',
                    line = list(color = '#1B9E77')) %>%
                layout(yaxis = list(title = "Total Log10 SARS CoV-2 Copies", 
                                 showline = TRUE,
                                 automargin = TRUE)) %>%
                layout(xaxis = list(title = "Date")) %>%
                layout(title = "WRF A") %>%
  plotly::add_markers(x = ~date, y = ~log_total_copies_both,
                      data = wrfa_both,
                       hoverinfo = "text",
                       showlegend = FALSE,
                       text = ~paste('</br> Date: ', date, 
                                     '</br> Actual Log Copies: ', round(log_total_copies_both, digits = 2)),
                       marker = list(color = '#1B9E77', size = 6, opacity = 0.65))

p_wrf_a
save(p_wrf_a, file = "./plotly_objs/p_wrf_a.rda")
#**************************************WRF B PLOT**********************************************
#add trendlines 
#extract data from geom_smooth
#both extract
# *********************************span 0.6***********************************
#*****************Must always update the n = TOTAL NUMBER OF DAYS*************************
extract_bothb <- ggplot(wrfb_both, aes(x = date, y = log_total_copies_both)) + 
  stat_smooth(aes(outfit=fit_bothb<<-..y..), method = "loess", color = '#D95F02', 
              span = 0.3, n = 457)
## Warning: Ignoring unknown aesthetics: outfit
#look at the fits to align dates and total observations
#both
extract_bothb
## `geom_smooth()` using formula 'y ~ x'

fit_bothb
##   [1] 12.51250 12.51353 12.51454 12.51553 12.51651 12.51748 12.51845 12.51943
##   [9] 12.52043 12.52144 12.52247 12.52353 12.52463 12.52578 12.52696 12.52821
##  [17] 12.52951 12.53088 12.53231 12.53383 12.53543 12.53711 12.53887 12.54067
##  [25] 12.54252 12.54439 12.54630 12.54823 12.55017 12.55213 12.55408 12.55604
##  [33] 12.55800 12.55994 12.56186 12.56376 12.56562 12.56746 12.56925 12.57099
##  [41] 12.57268 12.57431 12.57588 12.57740 12.57888 12.58033 12.58175 12.58314
##  [49] 12.58451 12.58585 12.58718 12.58849 12.58978 12.59106 12.59232 12.59358
##  [57] 12.59484 12.59609 12.59733 12.59858 12.59983 12.60109 12.60236 12.60363
##  [65] 12.60492 12.60622 12.60754 12.60888 12.61024 12.61163 12.61304 12.61462
##  [73] 12.61651 12.61868 12.62108 12.62370 12.62651 12.62946 12.63255 12.63573
##  [81] 12.63897 12.64225 12.64554 12.64880 12.65201 12.65514 12.65815 12.66102
##  [89] 12.66373 12.66623 12.66850 12.67051 12.67127 12.66995 12.66670 12.66168
##  [97] 12.65506 12.64700 12.63767 12.62722 12.61582 12.60363 12.59081 12.57754
## [105] 12.56397 12.55026 12.53658 12.52309 12.50995 12.49733 12.48539 12.47429
## [113] 12.46419 12.45527 12.44768 12.44158 12.43714 12.43451 12.43388 12.43539
## [121] 12.43860 12.44290 12.44818 12.45434 12.46128 12.46889 12.47708 12.48574
## [129] 12.49478 12.50407 12.51354 12.52306 12.53255 12.54394 12.55888 12.57686
## [137] 12.59735 12.61983 12.64377 12.66864 12.69393 12.71911 12.74365 12.76704
## [145] 12.78874 12.80823 12.82499 12.84278 12.86488 12.89012 12.91730 12.94524
## [153] 12.97275 12.99864 13.02172 13.04080 13.05890 13.07952 13.10218 13.12642
## [161] 13.15175 13.17770 13.20379 13.22955 13.25450 13.27817 13.30009 13.31977
## [169] 13.33674 13.35054 13.36322 13.37703 13.39174 13.40712 13.42295 13.43898
## [177] 13.45500 13.47077 13.48606 13.50065 13.51431 13.52680 13.53790 13.54737
## [185] 13.55500 13.56054 13.56377 13.56446 13.56239 13.55779 13.55116 13.54268
## [193] 13.53251 13.52080 13.50771 13.49342 13.47807 13.46184 13.44488 13.42737
## [201] 13.40945 13.39129 13.37305 13.35490 13.33699 13.31741 13.29448 13.26875
## [209] 13.24075 13.21102 13.18011 13.14855 13.11688 13.08564 13.05537 13.02661
## [217] 12.99990 12.97159 12.93820 12.90057 12.85951 12.81584 12.77040 12.72399
## [225] 12.67746 12.63161 12.58727 12.54527 12.50642 12.47155 12.44149 12.41095
## [233] 12.37521 12.33591 12.29473 12.25330 12.21329 12.17635 12.14413 12.11828
## [241] 12.09647 12.07528 12.05475 12.03492 12.01583 11.99753 11.98006 11.96345
## [249] 11.94775 11.93301 11.91926 11.90654 11.89491 11.88439 11.87503 11.86687
## [257] 11.85996 11.85434 11.85004 11.84711 11.84559 11.84781 11.85549 11.86781
## [265] 11.88394 11.90303 11.92425 11.94678 11.96977 11.99238 12.01380 12.03318
## [273] 12.04968 12.06248 12.07073 12.07903 12.09187 12.10830 12.12736 12.14808
## [281] 12.16951 12.19068 12.21064 12.22843 12.24308 12.25365 12.25916 12.26244
## [289] 12.26674 12.27175 12.27716 12.28266 12.28792 12.29265 12.29652 12.29923
## [297] 12.30046 12.29990 12.29723 12.29215 12.28435 12.27166 12.25311 12.23011
## [305] 12.20409 12.17645 12.14863 12.12203 12.09809 12.07822 12.05795 12.03245
## [313] 12.00258 11.96923 11.93327 11.89559 11.85704 11.81852 11.78090 11.74505
## [321] 11.71186 11.68219 11.65692 11.63694 11.61901 11.59956 11.57896 11.55760
## [329] 11.53584 11.51406 11.49263 11.47191 11.45229 11.43414 11.41783 11.40372
## [337] 11.39221 11.38160 11.37020 11.35832 11.34627 11.33438 11.32295 11.31231
## [345] 11.30276 11.29463 11.28824 11.28389 11.28191 11.28260 11.28470 11.28691
## [353] 11.28951 11.29282 11.29711 11.30268 11.30984 11.31886 11.33005 11.34434
## [361] 11.36226 11.38346 11.40761 11.43438 11.46342 11.49441 11.52700 11.56086
## [369] 11.59565 11.63103 11.66668 11.70225 11.73740 11.77181 11.80512 11.83702
## [377] 11.86716 11.89520 11.92455 11.95804 11.99452 12.03278 12.07167 12.11000
## [385] 12.14661 12.18030 12.20990 12.23875 12.27059 12.30487 12.34102 12.37850
## [393] 12.41676 12.45523 12.49337 12.53063 12.56644 12.60025 12.63152 12.65968
## [401] 12.68418 12.70639 12.72803 12.74904 12.76938 12.78903 12.80792 12.82603
## [409] 12.84331 12.85971 12.87520 12.88974 12.90328 12.91598 12.92799 12.93930
## [417] 12.94989 12.95974 12.96883 12.97715 12.98467 12.99138 12.99726 13.00228
## [425] 13.00643 13.00970 13.01206 13.01359 13.01439 13.01443 13.01369 13.01215
## [433] 13.00979 13.00660 13.00254 12.99761 12.99183 12.98524 12.97784 12.96964
## [441] 12.96064 12.95083 12.94022 12.92881 12.91659 12.90358 12.88977 12.87516
## [449] 12.85975 12.84355 12.82655 12.80876 12.79018 12.77080 12.75063 12.72967
## [457] 12.70792
#assign fits to a vector
both_trendb <- fit_bothb

#extract y min and max for each
limits_bothb <- ggplot_build(extract_bothb)$data
## `geom_smooth()` using formula 'y ~ x'
limits_bothb <- as.data.frame(limits_bothb)
both_yminb <- limits_bothb$ymin
both_ymaxb <- limits_bothb$ymax

#reassign dataframes (just to be safe)
work_bothb <- wrfb_both

#fill in missing dates to smooth fits
work_bothb <- work_bothb %>% complete(date = seq(min(date), max(date), by = "1 day"))
date_vec_bothb <- work_bothb$date

#create a new smooth dataframe to layer
smooth_frame_bothb <- data.frame(date_vec_bothb, both_trendb, both_yminb, both_ymaxb)
#WRF B
#plot smooth frames
p_wrf_b <- plotly::plot_ly() %>%
  plotly::add_lines(x = ~date_vec_bothb, y = ~both_trendb,
                    data = smooth_frame_bothb,
                    hoverinfo = "text",
                    text = ~paste('</br> Date: ', date_vec_bothb,
                                  '</br> Median Log Copies: ', round(both_trendb, digits = 2)),
                    line = list(color = '#D95F02', size = 8, opacity = 0.65),
                    showlegend = FALSE) %>%
     layout(xaxis = list(range = c(mindate - 7, maxdate + 7))) %>% #buffer here
plotly::add_ribbons(x ~date_vec_bothb, ymin = ~both_yminb, ymax = ~both_ymaxb,
                    showlegend = FALSE,
                    opacity = 0.25,
                    hoverinfo = "text",
                    text = ~paste('</br> Date: ', date_vec_bothb, #leaving in case we want to change
                                  '</br> Max Log Copies: ', round(both_ymaxb, digits = 2),
                                  '</br> Min Log Copies: ', round(both_yminb, digits = 2)),
                    name = "",
                    fillcolor = '#D95F02',
                    line = list(color = '#D95F02')) %>%
                layout(yaxis = list(title = "Total Log10 SARS CoV-2 Copies", 
                                 showline = TRUE,
                                 automargin = TRUE)) %>%
                layout(xaxis = list(title = "Date")) %>%
                layout(title = "WRF B") %>%
  plotly::add_markers(x = ~date, y = ~log_total_copies_both,
                      data = wrfb_both,
                       hoverinfo = "text",
                       showlegend = FALSE,
                       text = ~paste('</br> Date: ', date, 
                                     '</br> Actual Log Copies: ', round(log_total_copies_both, digits = 2)),
                       marker = list(color = '#D95F02', size = 6, opacity = 0.65))

p_wrf_b
save(p_wrf_b, file = "./plotly_objs/p_wrf_b.rda")

#**************************************WRF C PLOT********************************************** #add trendlines #extract data from geom_smooth # *********************************span 0.6*********************************** #*****************Must always update the n = TOTAL NUMBER OF DAYS*************************

extract_bothc <- ggplot(wrfc_both, aes(x = date, y = log_total_copies_both)) + 
  stat_smooth(aes(outfit=fit_bothc<<-..y..), method = "loess", color = '#E7298A', 
              span = 0.3, n = 457)
## Warning: Ignoring unknown aesthetics: outfit
#look at the fits to align dates and total observations
#both
extract_bothc
## `geom_smooth()` using formula 'y ~ x'

fit_bothc
##   [1] 11.54661 11.56730 11.58773 11.60788 11.62774 11.64729 11.66653 11.68544
##   [9] 11.70401 11.72222 11.74007 11.75755 11.77463 11.79132 11.80759 11.82344
##  [17] 11.83885 11.85381 11.86830 11.88233 11.89586 11.90890 11.92143 11.93344
##  [25] 11.94491 11.95584 11.96621 11.97600 11.98521 11.99381 12.00176 12.00911
##  [33] 12.01589 12.02213 12.02785 12.03309 12.03788 12.04225 12.04623 12.04984
##  [41] 12.05313 12.05612 12.05883 12.06131 12.06358 12.06568 12.06762 12.06945
##  [49] 12.07119 12.07287 12.07390 12.07371 12.07236 12.06994 12.06651 12.06215
##  [57] 12.05694 12.05094 12.04422 12.03687 12.02896 12.02055 12.01173 12.00255
##  [65] 11.99311 11.98347 11.97370 11.96387 11.95407 11.94436 11.93482 11.92552
##  [73] 11.91653 11.90793 11.89978 11.89217 11.88517 11.87884 11.87144 11.86136
##  [81] 11.84889 11.83432 11.81796 11.80010 11.78103 11.76104 11.74043 11.71950
##  [89] 11.69854 11.67785 11.65772 11.63844 11.62031 11.60362 11.58867 11.57576
##  [97] 11.56517 11.55721 11.55217 11.54834 11.54389 11.53893 11.53356 11.52787
## [105] 11.52197 11.51595 11.50991 11.50396 11.49820 11.49272 11.48763 11.48303
## [113] 11.47901 11.47568 11.47313 11.47147 11.47080 11.47121 11.47282 11.47571
## [121] 11.47998 11.48575 11.49310 11.50214 11.51297 11.52568 11.54313 11.56700
## [129] 11.59566 11.62749 11.66087 11.69418 11.72578 11.75405 11.77737 11.80028
## [137] 11.82794 11.85961 11.89456 11.93206 11.97138 12.01179 12.05256 12.09297
## [145] 12.13228 12.16976 12.20468 12.23632 12.26393 12.29055 12.31934 12.34981
## [153] 12.38145 12.41378 12.44629 12.47850 12.50991 12.54001 12.56832 12.59434
## [161] 12.61757 12.64002 12.66384 12.68869 12.71427 12.74025 12.76633 12.79217
## [169] 12.81748 12.84192 12.86519 12.88697 12.90693 12.92477 12.94016 12.95411
## [177] 12.96779 12.98113 12.99406 13.00654 13.01848 13.02983 13.04052 13.05049
## [185] 13.05968 13.06803 13.07546 13.08191 13.08733 13.09165 13.09480 13.09672
## [193] 13.09734 13.09661 13.09446 13.09083 13.08594 13.08006 13.07318 13.06529
## [201] 13.05638 13.04644 13.03547 13.02345 13.01037 12.99623 12.98102 12.96472
## [209] 12.94733 12.92884 12.90889 12.88730 12.86427 12.84000 12.81470 12.78858
## [217] 12.76186 12.73473 12.70741 12.67659 12.63949 12.59713 12.55051 12.50064
## [225] 12.44854 12.39520 12.34163 12.28885 12.23786 12.18967 12.14529 12.10573
## [233] 12.07199 12.03734 11.99541 11.94777 11.89598 11.84162 11.78624 11.73141
## [241] 11.67870 11.62968 11.58591 11.54896 11.52039 11.49767 11.47715 11.45873
## [249] 11.44227 11.42769 11.41485 11.40366 11.39399 11.38575 11.37881 11.37306
## [257] 11.36840 11.36470 11.36186 11.35977 11.35832 11.35739 11.35687 11.35665
## [265] 11.35661 11.35666 11.36112 11.37356 11.39271 11.41725 11.44593 11.47744
## [273] 11.51050 11.54382 11.57612 11.60611 11.63251 11.65403 11.66938 11.67728
## [281] 11.68518 11.70000 11.71958 11.74179 11.76448 11.78551 11.80275 11.81403
## [289] 11.81724 11.81397 11.80751 11.79820 11.78638 11.77241 11.75662 11.73937
## [297] 11.72099 11.70184 11.68226 11.66259 11.64319 11.62440 11.60656 11.58311
## [305] 11.54868 11.50546 11.45561 11.40130 11.34471 11.28801 11.23338 11.18297
## [313] 11.13898 11.10356 11.07888 11.06004 11.04079 11.02136 11.00196 10.98281
## [321] 10.96413 10.94612 10.92902 10.91303 10.89837 10.88526 10.87391 10.86455
## [329] 10.85738 10.85283 10.85086 10.85096 10.85266 10.85546 10.85888 10.86243
## [337] 10.86561 10.86794 10.87058 10.87495 10.88095 10.88849 10.89745 10.90774
## [345] 10.91926 10.93191 10.94559 10.96019 10.97562 10.99178 11.00857 11.02588
## [353] 11.04362 11.06168 11.07997 11.09838 11.11681 11.13517 11.15336 11.17234
## [361] 11.19304 11.21531 11.23895 11.26382 11.28973 11.31651 11.34400 11.37202
## [369] 11.40041 11.42899 11.45759 11.48605 11.51418 11.54183 11.56883 11.59499
## [377] 11.62015 11.64415 11.66713 11.68949 11.71140 11.73303 11.75455 11.77614
## [385] 11.79796 11.82018 11.84299 11.86794 11.89609 11.92691 11.95984 11.99433
## [393] 12.02983 12.06578 12.10163 12.13683 12.17084 12.20310 12.23305 12.26015
## [401] 12.28384 12.30597 12.32856 12.35137 12.37412 12.39657 12.41845 12.43950
## [409] 12.45946 12.47808 12.49509 12.51023 12.52325 12.53482 12.54578 12.55608
## [417] 12.56568 12.57455 12.58264 12.58992 12.59636 12.60190 12.60652 12.61017
## [425] 12.61282 12.61442 12.61495 12.61458 12.61351 12.61169 12.60906 12.60559
## [433] 12.60122 12.59591 12.58960 12.58225 12.57391 12.56466 12.55449 12.54340
## [441] 12.53141 12.51850 12.50468 12.48996 12.47433 12.45779 12.44034 12.42199
## [449] 12.40274 12.38258 12.36152 12.33956 12.31670 12.29294 12.26828 12.24273
## [457] 12.21628
#assign fits to a vector
both_trendc <- fit_bothc

#extract y min and max for each
limits_bothc <- ggplot_build(extract_bothc)$data
## `geom_smooth()` using formula 'y ~ x'
limits_bothc <- as.data.frame(limits_bothc)
both_yminc <- limits_bothc$ymin
both_ymaxc <- limits_bothc$ymax

#reassign dataframes (just to be safe)
work_bothc <- wrfc_both

#fill in missing dates to smooth fits
work_bothc <- work_bothc %>% complete(date = seq(min(date), max(date), by = "1 day"))
date_vec_bothc <- work_bothc$date

#create a new smooth dataframe to layer
smooth_frame_bothc <- data.frame(date_vec_bothc, both_trendc, both_yminc, both_ymaxc)
#WRF C
#plot smooth frames
p_wrf_c <- plotly::plot_ly() %>%
  plotly::add_lines(x = ~date_vec_bothc, y = ~both_trendc,
                    data = smooth_frame_bothc,
                    hoverinfo = "text",
                    text = ~paste('</br> Date: ', date_vec_bothc,
                                  '</br> Median Log Copies: ', round(both_trendc, digits = 2)),
                    line = list(color = '#E7298A', size = 8, opacity = 0.65),
                    showlegend = FALSE) %>%
     layout(xaxis = list(range = c(mindate - 7, maxdate + 7))) %>% #buffer here
plotly::add_ribbons(x ~date_vec_bothc, ymin = ~both_yminc, ymax = ~both_ymaxc,
                    showlegend = FALSE,
                    opacity = 0.25,
                    hoverinfo = "text",
                    text = ~paste('</br> Date: ', date_vec_bothc, #leaving in case we want to change
                                  '</br> Max Log Copies: ', round(both_ymaxc, digits = 2),
                                  '</br> Min Log Copies: ', round(both_yminc, digits = 2)),
                    name = "",
                    fillcolor = '#E7298A',
                    line = list(color = '#E7298A')) %>%
                layout(yaxis = list(title = "Total Log10 SARS CoV-2 Copies", 
                                 showline = TRUE,
                                 automargin = TRUE)) %>%
                layout(xaxis = list(title = "Date")) %>%
                layout(title = "WRF C") %>%
  plotly::add_markers(x = ~date, y = ~log_total_copies_both,
                      data = wrfc_both,
                       hoverinfo = "text",
                       showlegend = FALSE,
                       text = ~paste('</br> Date: ', date, 
                                     '</br> Actual Log Copies: ', round(log_total_copies_both, digits = 2)),
                       marker = list(color = '#E7298A', size = 6, opacity = 0.65))

p_wrf_c
save(p_wrf_c, file = "./plotly_objs/p_wrf_c.rda")
save(wrfa_both, file = "./plotly_objs/wrfa_both.rda")
save(wrfb_both, file = "./plotly_objs/wrfb_both.rda")
save(wrfc_both, file = "./plotly_objs/wrfc_both.rda")
save(date_vec_botha, file = "./plotly_objs/date_vec_botha.rda")
save(date_vec_bothb, file = "./plotly_objs/date_vec_bothb.rda")
save(date_vec_bothc, file = "./plotly_objs/date_vec_bothc.rda")
save(both_ymina, file = "./plotly_objs/both_ymina.rda")
save(both_ymaxa, file = "./plotly_objs/both_ymaxa.rda")

save(both_yminb, file = "./plotly_objs/both_yminb.rda")
save(both_ymaxb, file = "./plotly_objs/both_ymaxb.rda")

save(both_yminc, file = "./plotly_objs/both_yminc.rda")
save(both_ymaxc, file = "./plotly_objs/both_ymaxc.rda")